%% import data from .out files
% %
% label='1114';
% [all_wind,all_rb,fir_wind,fir_rb]=JOBdata('1114');
% 1  2   3  4   5  6   7   8   9  10 11 12  13    14 15 16  17  18 19
% m1,m2,Lx,ecc0,a,tb0,mt1,mt2,kw,kw2,mx,mx2,eccx,tbx,f,nnn,uuu,m_tr,i
% m1,m2,Lx,ecc0,a,tb0,mt1,mt2,kw,kw2,mx,mx2,eccx,tbx,f,m_transferrate,b_iso,m_tr,i
% 20 \delta T



%% initial parameters
global tmax ge AllmassNUM M1MAX M1MIN SFR qNUM aNUM ...
    minLx maxLx alpha CHAbeamingFlag eddfac

alpha=2.7;
tmax=1d4*1d6; % in units of year
AllmassNUM=1d3;

% qNUM=10;
aNUM=100;
M1MAX=150;
M1MIN=2;
SFR=GetSFR(50,alpha,0.535,0,M1MIN,5);
% SFR=10.6984;
yearsc=3.1557d+07;
Msun=1.9891d30;
minLx=39;
maxLx=41;
% xi=Mdot/(sum(sum(data(:,1:2)))/tmax)/tmax*1d6;
% warning!! 
CHAbeamingFlag=0;
eddfac=10000;
ge=30;
%% delete some useless data

[all_rb_1,all_wind_1]=DELallrbwind(all_rb,all_wind);
% allinall=[all_rb_1;all_wind_1];
%% plot wiktorowicz 2017 figure 2
tm=1d4;
all_rb_1(all_rb_1(:,3)>1d42,:)=[];
BHsum=plot_fir_gedian_lx(all_rb_1(all_rb_1(:,10)==14,:),1d4);
NSsum=plot_fir_gedian_lx(all_rb_1(all_rb_1(:,10)==13,:),1d4);
clear opt
plot(1:tm,NSsum./BHsum)
opt.XLabel='T / Myr';
opt.YLabel='ratio';
opt.Title='constant';
opt.Legend={'E(N_{NS}) / E(N_{BH}'};
opt.YScale='linear';
opt.XScale='log';
opt.XLim=[1 tm];
opt.YLim=[0 1.3];
opt.LegendBox='on';
setPlotProp(opt);
hold off









BHsum=plot_fir_gedian(all_rb_1(all_rb_1(:,10)==14,:),100);

NSsum=plot_fir_gedian(all_rb_1(all_rb_1(:,10)==13,:),100);
plot(1:tm,BHsum)
hold on
plot(1:tm,NSsum)
clear opt
opt.XLabel='T / Myr';
opt.YLabel='Number of ULXs';
opt.Title='Burst';
opt.Legend={'BH','NS'};
opt.YScale='log';
opt.XScale='log';
opt.XLim=[1 tm];
opt.YLim=[1d2 max(BHsum)];
opt.LegendBox='on';
setPlotProp(opt);
hold off

tm=1d4;
BHsum=plot_fir_gedian(all_rb_1(all_rb_1(:,10)==14,:),1d4);

NSsum=plot_fir_gedian(all_rb_1(all_rb_1(:,10)==13,:),1d4);
plot(1:tm,BHsum)
hold on
plot(1:tm,NSsum)
opt.XLabel='T / Myr';
opt.YLabel='Number of ULXs';
opt.Title='constant';
opt.Legend={'BH','NS'};
opt.YScale='log';
opt.XScale='log';
opt.XLim=[1 tm];
opt.YLim=[1d2 max(BHsum)];
opt.LegendBox='on';
setPlotProp(opt);
hold off

